#pragma once
#include <cuda_runtime.h>
#include <Utilities/glm_ext.h>

namespace PhysLeo {

/**
 * FEM Neohookean material.
 * \f$ I_1 = tr(F^TF), J = det(F)  \f$,
 * \f$ W(E) = \mu(\frac{I_1-2log(J)-3}{2})+\frac{\lambda}{2}log^2(J) \f$
 */
template<typename T>
class FemNeohookean {
public:
    /**
    * return the energy density according to deformation gradient
    * @param[in] F  a glm 3x3 matrix type represents the deformation gradient.
    * @param[in] ptr_material  a pointer to vertex specific material data
    * @return a T value, the energy density.
    */
    __host__ __device__ static T energyDensity(glm::tmat3x3<T> F, T* ptr_material)
    {
        auto lambda = ptr_material[0];
        auto mu = ptr_material[1];

        glm::tmat3x3<T> FFT = F * glm::transpose(F);
        T I1 = glm::trace(FFT);
        T J = glm::determinant(F);

        //W = mu/2*(I1-3)-mu*log(J)+lambda/2*log(J)^2
        return 0.5f*(mu*(I1 - 3)) - mu * log(J) + 0.5f*lambda*log(J)*log(J);
    }

    /**
    * return the first piola stress according to deformation gradient.
    * @param[in] F  a glm 3x3 matrix type represents the deformation gradient.
    * @param[in] ptr_material  a pointer to vertex specific material data
    * @return a glm 3x3 matrix type represents the first piola stress.
    */
    __host__ __device__ static glm::tmat3x3<T> firstPiolaStress(glm::tmat3x3<T> F, T* ptr_material)
    {
        auto lambda = ptr_material[0];
        auto mu = ptr_material[1];

        glm::tmat3x3<T> F_T_I = glm::inverse(glm::transpose(F));

        //P = mu*(F-F^(-T))+lambda*log(J)F^(-T)
        T log_detF = log(glm::determinant(F));
        return mu*(F - F_T_I) + lambda* log_detF *F_T_I;
    }
};

}